{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Project 2: 全景图像拼接\n",
    "\n",
    "启动Jupyter并打开`Project2.ipynb`。\n",
    "\n",
    "在本项目中，你将实现一个全景图像拼接器。这将建立在项目 1 中开发的兴趣点和描述符的基础上。您将首先通过 RANSAC 进行几何过滤，然后估计成对旋转并将其链接在一起以对齐全景图。当您有一个基本的缝合器工作时，使用更好的对齐、混合或其他新功能来改进它，并记录您的发现。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.linalg\n",
    "import os.path\n",
    "from time import time\n",
    "import types\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import im_util\n",
    "import interest_point\n",
    "import ransac\n",
    "import geometry\n",
    "import render\n",
    "import panorama\n",
    "\n",
    "%matplotlib inline\n",
    "# edit this line to change the figure size\n",
    "plt.rcParams['figure.figsize'] = (16.0, 10.0)\n",
    "# force auto-reload of import modules before running code \n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 图像变形测试\n",
    "\n",
    "下面的代码使用一个3x3变换矩阵对图像进行变形。使用矩阵P进行实验，测试课堂上描述的一些不同的二维变换，例如，相似性变换、仿射变换和投影变换。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read image\n",
    "image_dir='data/test'\n",
    "im_filename1=image_dir+'/100-0038_img.jpg'\n",
    "im=im_util.image_open(im_filename1)\n",
    "im_rows,im_cols,_=im.shape\n",
    "\n",
    "# set transformation matrix\n",
    "P=[[1, 0.2, -64],\n",
    "  [ 0, 1.1, -120],\n",
    "  [ 0, 5.2e-4, 0.83]]\n",
    "\n",
    "# warp coordinates\n",
    "r0,r1=-im_rows/2, im_rows*3/2\n",
    "c0,c1=-im_cols/2, im_cols*3/2\n",
    "warp_rows, warp_cols=im_rows, im_cols\n",
    "\n",
    "coords=im_util.coordinate_image(warp_rows,warp_cols,r0,r1,c0,c1)\n",
    "coords_t=im_util.transform_coordinates(coords, P)\n",
    "\n",
    "# visualise result\n",
    "warp_im1=im_util.warp_image(im,coords)\n",
    "warp_im2=im_util.warp_image(im,coords_t)\n",
    "alpha=im_util.warp_image(np.ones((im_rows,im_cols,1)),coords_t)\n",
    "result_im=warp_im2*alpha + 0.5*warp_im1*(1-alpha)\n",
    "\n",
    "ax1=plt.subplot(1,1,1)\n",
    "plt.imshow(result_im)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 兴趣点测试\n",
    "\n",
    "我们将使用在项目1(\"Project 1\")中实现的兴趣点和描述符。如果你在让这些工作起来时遇到了困难，或者想要测试另一种匹配算法，你可以尝试使用 vlfeat/sift，见下文。注意你需要安装 cyvlfeat，例如，使用 conda 安装：`conda install -c menpo cyvlfeat`。\n",
    "\n",
    "```python\n",
    "from cyvlfeat import sift \n",
    "frames, desc = sift.sift(img, compute_descriptor=True, n_levels=1)\n",
    "ip = (frames.T)[0:2, :]\n",
    "desc = desc.astype(np.float)\n",
    "```\n",
    "\n",
    "运行下面的两个代码块来检查你的兴趣点和描述符是否工作正常。为了让后续步骤顺利进行，你应该争取获得大约100-1000个兴趣点。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Read a pair of input images and convert to grey\n",
    "\"\"\"\n",
    "image_dir='data/test'\n",
    "#im_filename1=image_dir+'/100-0023_img.jpg'\n",
    "#im_filename2=image_dir+'/100-0024_img.jpg'\n",
    "im_filename1=image_dir+'/100-0038_img.jpg'\n",
    "im_filename2=image_dir+'/100-0039_img.jpg'\n",
    "\n",
    "im1 = im_util.image_open(im_filename1)\n",
    "im2 = im_util.image_open(im_filename2)\n",
    "\n",
    "img1 = np.mean(im1, 2, keepdims=True)\n",
    "img2 = np.mean(im2, 2, keepdims=True)\n",
    "\n",
    "#optionally plot images\n",
    "#ax1,ax2=im_util.plot_two_images(im1, im2)\n",
    "\n",
    "\"\"\"\n",
    "Find interest points in the image pair\n",
    "\"\"\"\n",
    "print('[ find interest points ]')\n",
    "t0=time()\n",
    "ip_ex = interest_point.InterestPointExtractor()\n",
    "ip1 = ip_ex.find_interest_points(img1)\n",
    "print(' found '+str(ip1.shape[1])+' in image 1')\n",
    "ip2 = ip_ex.find_interest_points(img2)\n",
    "print(' found '+str(ip2.shape[1])+' in image 2')\n",
    "t1=time()\n",
    "print(' % .2f secs ' % (t1-t0))\n",
    "\n",
    "# optionally draw interest points\n",
    "#print('[ drawing interest points ]')\n",
    "#ax1,ax2=im_util.plot_two_images(im1,im2)\n",
    "#interest_point.draw_interest_points_ax(ip1, ax1)\n",
    "#interest_point.draw_interest_points_ax(ip2, ax2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Extract and match descriptors\n",
    "\"\"\"\n",
    "print('[ extract descriptors ]')\n",
    "t0=time()\n",
    "desc_ex = interest_point.DescriptorExtractor()\n",
    "desc1 = desc_ex.get_descriptors(img1, ip1)\n",
    "desc2 = desc_ex.get_descriptors(img2, ip2)\n",
    "t1=time()\n",
    "print(' % .2f secs' % (t1-t0))\n",
    "\n",
    "print('[ match descriptors ]')\n",
    "t0=time()\n",
    "match_idx = desc_ex.match_descriptors(desc1, desc2)\n",
    "t1=time()\n",
    "print(' % .2f secs' % (t1-t0))\n",
    "\n",
    "ipm=ip2[:,match_idx]\n",
    "\n",
    "print('[ drawing matches ]')\n",
    "t0=time()\n",
    "ax1,ax2=im_util.plot_two_images(im1,im2)\n",
    "interest_point.draw_matches_ax(ip1, ipm, ax1, ax2)\n",
    "t1=time()\n",
    "print(' % .2f secs' % (t1-t0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RANSAC 实现\n",
    "\n",
    "我们现在将使用 RANSAC 来找到一致的匹配。\n",
    "\n",
    "### 一致性测试 [3分]\n",
    "\n",
    "首先，我们将实现一个测试，以计算与相似性变换一致的匹配数量。下面的代码生成一个随机的相似性变换 S 和一组随机点 x。然后它变换这些点并添加噪声，并检查有多少点与真实的相似性变换 S 一致。\n",
    "\n",
    "打开 `ransac.py` 并实现函数 `consistent`。当运行下面的代码时，你应该发现有高比例（约80%或更多）的点与真实的相似性变换 S 一致。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Test RANSAC functions using synthetic data\n",
    "\"\"\"\n",
    "# make a random S matrix\n",
    "sd_pos=100\n",
    "sd_angle=np.pi\n",
    "theta=np.random.randn()*sd_angle\n",
    "tx=np.random.randn()*sd_pos\n",
    "ty=np.random.randn()*sd_pos\n",
    "ct=np.cos(theta)\n",
    "st=np.sin(theta)\n",
    "S=[[ct,st,tx],[-st,ct,ty],[0,0,1]]\n",
    "\n",
    "# generate random points\n",
    "num_points=100\n",
    "sd_points=20\n",
    "x = np.random.randn(2,num_points)*sd_points\n",
    "xh = geometry.hom(x)\n",
    "\n",
    "# transform points and add noise\n",
    "sd_noise=5\n",
    "yh = np.dot(S, xh)\n",
    "y = geometry.unhom(yh)\n",
    "yn = y + np.random.randn(2,num_points)*sd_noise\n",
    "\n",
    "print('[ Test of consistent ]')\n",
    "rn = ransac.RANSAC()\n",
    "inliers0=rn.consistent(S,x,yn)\n",
    "num_consistent=np.sum(inliers0)\n",
    "print(' number of points consistent with true S = '+str(num_consistent))\n",
    "if (num_consistent > 0.75*num_points):\n",
    "    print(' consistency check is working!')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 相似性求解器 [4分]\n",
    "\n",
    "现在选择2个点对应的样本，并计算与该对应的相似性变换。在 `ransac.py` 中实现 `compute_similarity` 函数，并运行下面的代码来计算内点的数量。尝试改变样本的索引，观察内点数量的变化。是否存在退化情况？如何检测这些情况？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('[ Test of compute_similarity ]')\n",
    "sample=[0,1]\n",
    "S1=rn.compute_similarity(x[:,sample],yn[:,sample])\n",
    "inliers1=rn.consistent(S1,x,yn)\n",
    "num_consistent=np.sum(inliers1)\n",
    "print(' number of points consistent with sample S = '+str(num_consistent))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RANSAC 循环 [3分]\n",
    "\n",
    "最后，通过完成 `ransac.py` 中的 `ransac_similarity` 函数来完成 RANSAC 的实现。完成后，您应该会发现大多数点都被标记为一致。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('[ Test of ransac_similarity ]')\n",
    "S2, inliers2=rn.ransac_similarity(x, yn)\n",
    "num_consistent=np.sum(inliers2)\n",
    "print(' number of points consistent with ransac S = '+str(num_consistent))\n",
    "if (num_consistent > 0.75*num_points):\n",
    "    print(' ransac succeeded!')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在我们将离开合成测试数据，使用上面输入图像对获取的兴趣点匹配运行相同的代码。请查看下面的代码，并检查输出是否合理。您应该会看到一组良好的几何一致的匹配点。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Perform RANSAC on interest point matches\n",
    "\"\"\"\n",
    "print('[ do ransac ]')\n",
    "t0=time()\n",
    "rn = ransac.RANSAC()\n",
    "S, inliers = rn.ransac_similarity(ip1,ipm)\n",
    "t1=time()\n",
    "num_inliers_s = np.sum(inliers)\n",
    "print(' found '+str(num_inliers_s)+' matches')\n",
    "print(' % .2f secs' % (t1-t0))\n",
    "\n",
    "ip1c = ip1[:, inliers]\n",
    "ipmc = ipm[:, inliers]\n",
    "\n",
    "print('[ drawing matches ]')\n",
    "t0=time()\n",
    "ax1,ax2=im_util.plot_two_images(im1,im2)\n",
    "interest_point.draw_matches_ax(ip1c, ipmc, ax1, ax2)\n",
    "t1=time()\n",
    "print(' % .2f secs' % (t1-t0))\n",
    "\n",
    "# optionally plot descriptors for matched points\n",
    "#inlier_id=np.flatnonzero(inliers)\n",
    "#match_id=match_idx[inlier_id]\n",
    "#interest_point.plot_matching_descriptors(desc1,desc2,inlier_id,match_id,plt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 旋转估计 [3分]\n",
    "\n",
    "下一个任务是估计图像之间的真实旋转。为了做到这一点，我们将猜测输入图像的视场角，并使用一个闭合形式的算法来估计旋转。打开 `geometry.py` 并完成 `compute_rotation` 函数的实现。您应该会发现大量的匹配点与您的旋转一致，并且成对的变形图像应该看起来合理。尝试调整视场角参数。对于这些图像来说，最佳的视场角是多少？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Estimate rotation matrix by least squares\n",
    "\"\"\"\n",
    "print('[ estimate rotation ]')\n",
    "t0=time()\n",
    "# Note: assume field of view of 45 degrees\n",
    "fov_degrees=45\n",
    "print(' assuming fov='+str(fov_degrees))\n",
    "K1 = geometry.get_calibration(im1.shape, fov_degrees)\n",
    "K2 = geometry.get_calibration(im2.shape, fov_degrees)\n",
    "R,H = geometry.compute_rotation(ip1c, ipmc, K1, K2)\n",
    "\n",
    "num_inliers_r = np.sum(rn.consistent(H, ip1, ipm))\n",
    "print(' num consistent with rotation = '+str(num_inliers_r))\n",
    "if (num_inliers_r>0.9 * num_inliers_s):\n",
    "    print(' compute rotation succeeded!')\n",
    "t1=time()\n",
    "print(' % .2f secs' % (t1-t0))\n",
    "    \n",
    "print('[ test pairwise warp ]')\n",
    "t0=time()\n",
    "im1_w, im2_w = render.pairwise_warp(im1, im2, H)\n",
    "_= im_util.plot_two_images(0.5*(im1+im2_w), 0.5*(im2+im1_w))\n",
    "t1=time()\n",
    "print(' % .2f secs' % (t1-t0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "以下代码将在球坐标系中渲染对齐的图像。请检查图像是否对齐良好。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Render 2 images in spherical coordinates\n",
    "\"\"\"\n",
    "images=[im1,im2]\n",
    "P1=K1\n",
    "P2=np.dot(K2,R)\n",
    "P_matrices=[P1,P2]\n",
    "\n",
    "render_params={}\n",
    "render_params['render_width']=800\n",
    "render_params['theta_min']=-45\n",
    "render_params['theta_max']=45\n",
    "render_params['phi_min']=-30\n",
    "render_params['phi_max']=30\n",
    "\n",
    "print ('[ render aligned images ]')\n",
    "t0=time()\n",
    "pano_im=render.render_spherical(images, P_matrices, render_params)\n",
    "t1=time()\n",
    "print(' % .2f secs' % (t1-t0))\n",
    "\n",
    "plt.plot()\n",
    "plt.imshow(pano_im)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "让我们添加更多的图像！`panorama.py` 中的 `PanoramaStitcher` 类的方法 `match_images` 接受一组图像作为输入，并将兴趣点和匹配代码封装在其中。请查看这个函数，并使用下面的代码在一组图像上进行测试。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Read a set of input images\n",
    "\"\"\"\n",
    "\n",
    "print('[ read images ]')\n",
    "image_dir='data/test'\n",
    "im_filenames=os.listdir(image_dir)\n",
    "im_filenames=[image_dir+'/'+fname for fname in im_filenames]\n",
    "\n",
    "#im_filenames=[]\n",
    "#im_filenames.append(image_dir+'/100-0023_img.jpg')\n",
    "#im_filenames.append(image_dir+'/100-0024_img.jpg')\n",
    "#im_filenames.append(image_dir+'/100-0038_img.jpg')\n",
    "#im_filenames.append(image_dir+'/100-0039_img.jpg')\n",
    "\n",
    "images=[]\n",
    "for fname in im_filenames:\n",
    "  images.append(im_util.image_open(fname))\n",
    "\n",
    "\"\"\"\n",
    "Stitch images\n",
    "\"\"\"\n",
    "stitch_params={}\n",
    "stitch_params['fov_degrees']=45\n",
    "\n",
    "num_images = len(im_filenames)\n",
    "print(' stitching '+str(num_images)+' images')\n",
    "\n",
    "pano=panorama.PanoramaStitcher(images, stitch_params)\n",
    "pano.match_images()\n",
    "\n",
    "print(' num_matches=')\n",
    "print(pano.num_matches)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 全景对齐 [3分]\n",
    "\n",
    "现在编写代码，通过将成对的旋转矩阵连接在一起，为每个图像计算一个旋转矩阵（假设第一张图像是单位旋转）。这段代码应该放在 `panorama.py` 中的 `align_panorama` 函数中。\n",
    "\n",
    "现在，您可以使用 `render` 方法将所有图像在球坐标系中拼接起来，如下所示的代码所示。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pano.align_panorama()\n",
    "\n",
    "render_params={}\n",
    "render_params['render_width']=800\n",
    "render_params['theta_min']=-45\n",
    "render_params['theta_max']=45\n",
    "render_params['phi_min']=-30\n",
    "render_params['phi_max']=30\n",
    "\n",
    "pano_im = pano.render(render_params)\n",
    "\n",
    "plt.plot()\n",
    "plt.imshow(pano_im)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 测试和改进全景拼接器 [20%]\n",
    "\n",
    "您现在应该已经完成了一个基本的全景图拼接器的实现。尝试使用几组不同的图像集进行测试，并记录结果中的任何问题/伪影。如何改进结果？列出可能的改进措施，并考虑您可能想要添加的新功能。现在实现其中一些改进/新功能，并在下面的笔记本中记录您的工作。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### TODO your improvements to the panorama stitcher\n",
    "\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
